#####################################################################################################
#
#   TITLE: International Intervention and the Rule of Law after Civil War: Replication Files (Part 2)
#   AUTHOR: Robert A. Blair
#
#####################################################################################################

# Clear

      rm(list = ls())


# Set directory

      #setwd("YOUR DIRECTORY")
      setwd("/Users/robertblair/Dropbox (Personal)/Documents/Research/Quant/Paper/Submissions/IO Acceptance/Analysis Replication Files")


# Load packages

      library(arm)
      library(data.table)
      library(MASS)
      library(foreign)
      


# Replicate Figure 1

      df = fread("r2_res_freqtype_pref_mlogit_margins.txt", skip = 1, header = T, check.names=T)
      
      rows = df[c(3,5,7,9,11)]$V1
      cols = c("State", "Local", "UNMIL")
      
      r2_c_pref_all = df[c(3,5,7,9,11),
                             c('r2_c_pref_all','r2_c_pref_all.1','r2_c_pref_all.2')]
      r2_c_pref_all = as.matrix(sapply(r2_c_pref_all, as.numeric))
      rownames(r2_c_pref_all) = rows
      colnames(r2_c_pref_all) = cols
      
      r2_c_pref_all.se = df[c(4,6,8,10,12),
                                c('r2_c_pref_all','r2_c_pref_all.1','r2_c_pref_all.2')]
      r2_c_pref_all.se = as.matrix(sapply(r2_c_pref_all.se, as.numeric))
      rownames(r2_c_pref_all.se) = rows
      colnames(r2_c_pref_all.se) = cols
      
      pdf("Figure1.pdf", 10, 7)
      par(family="Times")
      coefplot(as.numeric(r2_c_pref_all[,2]), as.numeric(r2_c_pref_all.se[,2]), CI=2, 
               ylab="Change in predicted probability", main = "", mar = c(1,5.1,5.1,1), ylim=c(-.3, .3),
               xlim = c(1,5.5),vertical=FALSE, pch.pts = 15, var.las = 1, cex.pts = 1.5, cex.var =1,
               varnames = c('            weekly\n            patrols', '            monthly\n            patrols',
                            '            occasional\n            patrols', '        interventions',
                            '            public\n            works'))
      coefplot(as.numeric(r2_c_pref_all[,1]), as.numeric(r2_c_pref_all.se[,1]),
               vertical=FALSE, pch.pts = 16, cex.pts = 1.5, cex.var =1,
               add=TRUE, offset = 0.1)
      coefplot(as.numeric(r2_c_pref_all[,3]), as.numeric(r2_c_pref_all.se[,3]),
               vertical=FALSE, pch.pts = 17, cex.pts = 1.5, cex.var =1,
               add=TRUE, offset = 0.2)
      legend('bottom', 
             legend = c("rely on informal", "rely on formal", "rely on UNMIL"),
             bty = "o", ncol = 3, inset = 0.05, cex = 0.9, pt.cex = 1, pch = c(15, 16, 17))
      dev.off() 

      
      
# Replicate Figure 2

      df_anyres = fread("r2_res_anyres_pref_mlogit_margins.txt", skip = 1, header = T, check.names=T)
      cols = c("State", "Local", "UNMIL")
      
      r2_c_pref_all = df_anyres[c(3),c('r2_c_pref_all','r2_c_pref_all.1','r2_c_pref_all.2')]
      colnames(r2_c_pref_all) = cols
      r2_c_pref_all = as.matrix(r2_c_pref_all)
      
      r2_c_pref_all.se = df_anyres[c(4),c('r2_c_pref_all','r2_c_pref_all.1','r2_c_pref_all.2')]
      colnames(r2_c_pref_all.se) = cols
      r2_c_pref_all.se = as.matrix(r2_c_pref_all.se)
      
      pdf("Figure2.pdf", 9, 5)
      par(family="Times")
      coefplot(as.numeric(r2_c_pref_all[,2]), as.numeric(r2_c_pref_all.se[,2]), CI=2, 
               ylab="Change in predicted probability", main = "", mar = c(1,5.1,5.1,1), ylim=c(-.15, .15), xlim= c(0.85,1.45),
               vertical=FALSE, pch.pts = 15, var.las = 1, cex.pts = 1.5, cex.var =1,
               varnames = c('', pretty.xlabels=TRUE))
      coefplot(as.numeric(r2_c_pref_all[,1]), as.numeric(r2_c_pref_all.se[,1]),
               vertical=FALSE, pch.pts = 16, cex.pts = 1.5, cex.var =1,
               add=TRUE, offset = 0.15)
      coefplot(as.numeric(r2_c_pref_all[,3]), as.numeric(r2_c_pref_all.se[,3]),
               vertical=FALSE, pch.pts = 17, cex.pts = 1.5, cex.var =1,
               add=TRUE, offset = 0.3)
      legend('bottom', 
             legend = c("rely on informal", "rely on formal",   "rely on UNMIL"),
             bty = "o", ncol = 3, inset = 0.005, cex = 0.9, pt.cex = 1, pch = c(15, 16, 17))
      dev.off()
      

      
# Define functions to replicate Figure 3

    # The following is a slightly modified version of the Corstange LISTIT function. All the changes have been marked with comments

    listit <- function(Y, X, map=NULL, method="BFGS") {
      ## Returns a matrix of parameters with a row for each list item
      ## and columns for each covariate.  The parameters are placed in
      ## their appropriate cells, and the remaining cells are filled
      ## with zeros.
      place.B <- function(param, map) {
        k <- nrow(map)
        n.x <- ncol(map)
        n.param <- length(param)
        B <- matrix(NA, nrow=k,
                    ncol=n.x)
        i.param <- 1
        for (i in 1:k)
          for (j in 1:n.x) {
            if (map[i,j] == TRUE) {
              B[i,j] <- param[i.param]
              i.param <- i.param + 1
            }
            else
              B[i,j] <- 0
          }
        return(B)
      }
      
      neg.LL <- function(param, Y.t, Y.c, X.t, X.c, map) {
        k <- nrow(map)
        B <- place.B(param=param, map=map)
        phi <- apply((1 + exp(-X.t %*%
                                t(B[2:k,])))^(-1)
                     , 1, sum)
        p <- as.vector(((1 + exp(-X.t %*%
                                   B[1,]))^(-1) +
                          phi) / k)
        q <- (1 + exp(-X.c %*% t(B[2:k,])))^(-1)
        list.LL <- sum(log(choose(k, Y.t)) +
                         Y.t * log(p) +
                         (k - Y.t) * log(1 - p))
        off.LL <- sum(Y.c * log(q) +
                        (1 - Y.c) * log(1 - q))
        joint.LL <- list.LL + off.LL
        return(-joint.LL)
      }
      
      neg.grad <-function(param, Y.t, Y.c, X.t, X.c, map) {
        k <- nrow(map)                      #number of list items
        n.x <- ncol(map)                    #number of covariates
        n.grad <- sum(map)                  #number of gradients to return
        B <- place.B(param=param, map=map)  #rows: list item, cols: parameters
        phi <- apply((1 + exp(-X.t %*%
                                t(B[2:k,])))^(-1)
                     , 1, sum)
        p <- as.vector(((1 + exp(-X.t %*%
                                   B[1,]))^(-1) +
                          phi) / k)
        ## Calculate the derivative of p with respect to the parameters
        ## Loop through the list items
        ## For each list item, get a temporary matrix of its predictors (this.X)
        ## Use (i.grad) and (this.param) indices to make sure the results go into
        ##   the appropriate columns in the (d.p) matrix, which has a number of
        ##   rows equal to the number of treatment observations, and a number of
        ##   columns equal to the number of gradients to be calculated (i.e., the
        ##   number of parameters we're trying to estimate).
        i.grad <- 0                         #gradient index
        d.p <- matrix(0, nrow=nrow(Y.t), ncol=n.grad)
        for (i in 1:k) {
          this.X <- X.t[,1]
          this.param <- sum(map[i,])
          for (j in 2:n.x)
            if (map[i,j] == TRUE)
              this.X <- cbind(this.X, X.t[,j]) #only the Xs for this list item
          e.vec <-
            as.vector(exp(-X.t %*% B[i,]) /
                        (k * (1 + exp(-X.t %*% B[i,]))^(2)))
          d.p[,(1 + i.grad):
                (i.grad + this.param)] <-
            this.X * e.vec
          i.grad <- i.grad + this.param
        }
        ## Calculate the contributions to the gradient calculations from the
        ## off items (i.e., the non-sensitive list items).
        ## Loop through the off items (items 2 to k).
        ## For each item, get a temporary matrix of its predictors (this.X).
        ## Use (i.grad) and (this.param) indices to make sure the results go into
        ##   the appropriate column in the (off) matrix.  (Note that the (off)    
        ##   matrix leaves a number of initial columns as zeros, corresponding to
        ##   the sensitive betas which do not enter these calculations.)  The
        ##   number of rows is the number of control group observations, and the
        ##   number of columns is the number of gradients to be calculated (i.e.,
        ##   the number of parameters we're trying to estimate).
        i.grad <- sum(map[1,])              #Start after the sensitive betas
        off <- matrix(0, nrow=nrow(Y.c), ncol=n.grad)
        for (i in 2:k) {
          this.X <- X.c[,1]
          this.param <- sum(map[i,])
          for (j in 2:n.x)
            if (map[i,j] == TRUE)
              this.X <- cbind(this.X, X.c[,j])
          q <- as.vector((1 + exp(-X.c %*%
                                    B[i,]))^(-1))
          off[,(1 + i.grad):(i.grad + this.param)] <-
            (Y.c[,(i - 1)] - q) * this.X
          i.grad <- i.grad + this.param
        }
        ## Calculate all the relevant components and add them together.
        ## (LL.t) is the component related to the treatment group.
        ## (LL.off) is the component related to the off items.
        ## (LL.joint) is the joint components.
        LL.t <- as.vector(Y.t * p^(-1)) * d.p +
          as.vector((k - Y.t) * (1 - p)^(-1)) * (-d.p)
        LL.t <- apply(LL.t, 2, sum)
        LL.off <- apply(off, 2, sum)
        LL.joint <- LL.t + LL.off
        return(-LL.joint)
      }
      
      ## Removes observations that have missing values from the
      ## working data.
      na.rm <- function(Y=Y, X=X) {
        k <- ncol(Y)                          #number of list items
        x.names <- colnames(X)                #names of the covariates
        X <- cbind(1,X)                       #append a constant to X
        colnames(X) <- c("constant", x.names) #names the columns
        ## D is a data frame.  The first column will be a vector of
        ## dummies indicating whether these observations will be used
        ## or not depending on whether or not they have NAs among the
        ## covariates.  The second column marks whether the observation
        ## is in the treatment group (ie, it gets the list) or the
        ## control group (it gets the off items individually).  The
        ## Y values follow, then the X values (with an appended constant).
        D <- as.data.frame(cbind(NA,NA,Y,X))
        colnames(D) <- c("to.use", "control",
                         colnames(Y),
                         colnames(X))
        ## If one of the off items is NA, puts NA in control, else a number
        D[,"control"] <-
          apply(D[,(4:(k + 2))], 1, sum)
        ## If control is not NA, it's a control variable, else it's NA...
        D[,"control"] <-
          ifelse(!is.na(D[,"control"]),
                 TRUE, D[,"control"])
        ## ...But if there is a list count, it's not a control variable
        D[,"control"] <-
          ifelse(!is.na(D[,3]),
                 FALSE, D[,"control"])
        ## Puts a count of the number of covariate NAs into to.use
        D[,"to.use"] <-
          apply(is.na(D[,((k + 3):ncol(D))]),
                1, sum)
        ## Mark to.use if there are no covariate NAs on this obs
        D[,"to.use"] <-
          ifelse((D[,"to.use"] != 0),
                 FALSE, TRUE)
        ## Mark to.use as FALSE if control observations have NA Y values
        D[,"to.use"] <-
          ifelse(is.na(D[,"control"]),
                 FALSE, D[,"to.use"])
        ## Take the subset of usable observations, stripping off
        ## the to.use and control columns
        D <- subset(D[3:ncol(D)], (D[,"to.use"] == TRUE))
        return(D)                             #return the usable Y and X values
      }
      
      ## Checks for the map defining which covariates predict which list items.
      ## If no map is supplied by the user, it assumes that every covariate in X
      ## predicts every list item.  If a map is supplied by a user, it appends
      ## an initial column of "TRUE" values to indicate that the constant is
      ## a covariate for every list item.
      if (is.null(map))
        map <- matrix(TRUE, nrow=ncol(Y),
                      ncol=(ncol(X) + 1))
      else map <- cbind(TRUE, map)
      ## Processes the data.  First it removes observations with NA
      ## values, then places into separate objects control and treatment
      ## group observations on Y and X.
      k <- ncol(Y)
      YX <- na.rm(Y=Y, X=X)
      YX.c <- subset(YX, is.na(YX[,1]))
      YX.t <- subset(YX, !is.na(YX[,1]))
      Y.c <- as.matrix(YX.c[,2:k])
      X.c <- as.matrix(YX.c[(k + 1):ncol(YX.c)])
      Y.t <- as.matrix(YX.t[,1])
      X.t <- as.matrix(YX.t[(k + 1):ncol(YX.t)])
      ## The number of parameters are the sum total TRUE values
      ## in the parameter map
      n.param <- sum(map)
      
      result <- optim(par=rep(0, n.param),
                      fn=neg.LL, gr=neg.grad,
                      hessian=TRUE,
                      method=method,
                      Y.t=Y.t, Y.c=Y.c,
                      X.t=X.t, X.c=X.c,
                      map=map)
      
      neg.se <- FALSE
      for (i in 1:ncol(X.t))
        if (is.nan(sqrt(diag(solve(result$hessian)))[i])) {
          neg.se <- TRUE
          break
        }
      
      output <- list(b.star=(result$par)[1:ncol(X.t)],
                     var=solve(result$hessian),
                     se=sqrt(diag(solve(result$hessian)))[1:ncol(X.t)],
                     vcov=solve(result$hessian)[1:ncol(X.t),1:ncol(X.t)], # added line
                     neg.se=neg.se,
                     converged=(result$convergence == 0),
                     n.treatment=nrow(Y.t),
                     n.control=nrow(Y.c),
                     items=k,
                     logL=-result$value,
                     deviance=2 * result$value,  ## ")" was removed and "," added to the end of this line
                     ## the following three lines were added:
                     b.control=(result$par)[(ncol(X.t)+1):length(result$par)],
                     se.control=sqrt(diag(solve(result$hessian)))[(ncol(X.t)+1):length(result$par)],
                     vcov.control=solve(result$hessian)[(ncol(X.t)+1):length(result$par),(ncol(X.t)+1):length(result$par)])
      
      class(output) <- "listit"
      return(output)
    }

    logistic <- function(x) exp(x)/(1+exp(x))
    
    monte.carlo<-function(freq, freq.not, n.draws, draws){
      freq.list <- freq.not.list <- freq.diff <- rep(NA, n.draws)
      for (d in 1:n.draws) {
        par.g <- draws[d, ]
        pred.freq.list <- na.exclude(logistic(freq %*% par.g))
        pred.freq.not.list <- na.exclude(logistic(freq.not %*% par.g))
        freq.list[d] <- mean(pred.freq.list)
        freq.not.list[d] <- mean(pred.freq.not.list)
        freq.diff[d] <- freq.list[d] - freq.not.list[d]
      }
      freq.diff
    }

    find.coefs<-function(sensitive, control1, control2, control3, vcov, vcov1, vcov2, vcov3, n.draws=10000){
      
      draws <- mvrnorm(n = n.draws, mu = sensitive, Sigma = vcov) #sensitive item
      freq.diff <- monte.carlo(weekly, weekly.not, n.draws, draws) 
      weekly.diff.mean <- mean(freq.diff)
      weekly.diff.se <- sd(freq.diff)
      freq.diff <- monte.carlo(monthly, monthly.not, n.draws, draws) 
      monthly.diff.mean <- mean(freq.diff)
      monthly.diff.se <- sd(freq.diff)
      freq.diff <- monte.carlo(occasional, occasional.not, n.draws, draws) 
      occasional.diff.mean <- mean(freq.diff)
      occasional.diff.se <- sd(freq.diff)
      freq.diff <- monte.carlo(interventions, interventions.not, n.draws, draws) 
      interventions.diff.mean <- mean(freq.diff)
      interventions.diff.se <- sd(freq.diff)
      freq.diff <- monte.carlo(public.works, public.works.not, n.draws, draws) 
      public.works.diff.mean <- mean(freq.diff)
      public.works.diff.se <- sd(freq.diff)
      
      draws <- mvrnorm(n = n.draws, mu = control1, Sigma = vcov1) #Call police
      freq.diff <- monte.carlo(weekly, weekly.not, n.draws, draws) 
      weekly.diff.mean1 <- mean(freq.diff)
      weekly.diff.se1 <- sd(freq.diff)
      freq.diff <- monte.carlo(monthly, monthly.not, n.draws, draws) 
      monthly.diff.mean1 <- mean(freq.diff)
      monthly.diff.se1 <- sd(freq.diff)
      freq.diff <- monte.carlo(occasional, occasional.not, n.draws, draws) 
      occasional.diff.mean1 <- mean(freq.diff)
      occasional.diff.se1 <- sd(freq.diff)
      freq.diff <- monte.carlo(interventions, interventions.not, n.draws, draws) 
      interventions.diff.mean1 <- mean(freq.diff)
      interventions.diff.se1 <- sd(freq.diff)
      freq.diff <- monte.carlo(public.works, public.works.not, n.draws, draws) 
      public.works.diff.mean1 <- mean(freq.diff, na.rm=TRUE)
      public.works.diff.se1 <- sd(freq.diff, na.rm=TRUE)
      
      draws <- mvrnorm(n = n.draws, mu = control2, Sigma = vcov2) #Call UNMIL
      freq.diff <- monte.carlo(weekly, weekly.not, n.draws, draws) 
      weekly.diff.mean2 <- mean(freq.diff)
      weekly.diff.se2 <- sd(freq.diff)
      freq.diff <- monte.carlo(monthly, monthly.not, n.draws, draws) 
      monthly.diff.mean2 <- mean(freq.diff)
      monthly.diff.se2 <- sd(freq.diff)
      freq.diff <- monte.carlo(occasional, occasional.not, n.draws, draws) 
      occasional.diff.mean2 <- mean(freq.diff)
      occasional.diff.se2 <- sd(freq.diff)
      freq.diff <- monte.carlo(interventions, interventions.not, n.draws, draws) 
      interventions.diff.mean2 <- mean(freq.diff)
      interventions.diff.se2 <- sd(freq.diff)
      freq.diff <- monte.carlo(public.works, public.works.not, n.draws, draws) 
      public.works.diff.mean2 <- mean(freq.diff)
      public.works.diff.se2 <- sd(freq.diff)
      
      draws <- mvrnorm(n = n.draws, mu = control3, Sigma = vcov3) #Do nothing
      freq.diff <- monte.carlo(weekly, weekly.not, n.draws, draws) 
      weekly.diff.mean3 <- mean(freq.diff)
      weekly.diff.se3 <- sd(freq.diff)
      freq.diff <- monte.carlo(monthly, monthly.not, n.draws, draws) 
      monthly.diff.mean3 <- mean(freq.diff)
      monthly.diff.se3 <- sd(freq.diff)
      freq.diff <- monte.carlo(occasional, occasional.not, n.draws, draws) 
      occasional.diff.mean3 <- mean(freq.diff)
      occasional.diff.se3 <- sd(freq.diff)
      freq.diff <- monte.carlo(interventions, interventions.not, n.draws, draws) 
      interventions.diff.mean3 <- mean(freq.diff, na.rm=TRUE)
      interventions.diff.se3 <- sd(freq.diff, na.rm=TRUE)
      freq.diff <- monte.carlo(public.works, public.works.not, n.draws, draws) 
      public.works.diff.mean3 <- mean(freq.diff, na.rm=TRUE)
      public.works.diff.se3 <- sd(freq.diff, na.rm=TRUE)
  
      means <- data.frame(c(weekly.diff.mean, weekly.diff.mean1, weekly.diff.mean2, weekly.diff.mean3),
                          c(monthly.diff.mean, monthly.diff.mean1, monthly.diff.mean2, monthly.diff.mean3),
                          c(occasional.diff.mean, occasional.diff.mean1, occasional.diff.mean2, occasional.diff.mean3),
                          c(interventions.diff.mean, interventions.diff.mean1, interventions.diff.mean2, interventions.diff.mean3),
                          c(public.works.diff.mean, public.works.diff.mean1, public.works.diff.mean2, public.works.diff.mean3))
      rownames(means) <- cbind('Use trial by ordeal', 'Call police', 'Call UNMIL', 'Do nothing')
      colnames(means) <- c('weekly', 'monthly', 'occasional', 'interventions', 'public works')
      
      ses <- data.frame(c(weekly.diff.se, weekly.diff.se1, weekly.diff.se2, weekly.diff.se3),
                        c(monthly.diff.se, monthly.diff.se1, monthly.diff.se2, monthly.diff.se3),
                        c(occasional.diff.se, occasional.diff.se1, occasional.diff.se2, occasional.diff.se3),
                        c(interventions.diff.se, interventions.diff.se1, interventions.diff.se2, interventions.diff.se3),
                        c(public.works.diff.se, public.works.diff.se1, public.works.diff.se2, public.works.diff.se3))
      rownames(ses) <- cbind('Use trial by ordeal', 'Call police', 'Call UNMIL', 'Do nothing')
      colnames(ses) <- c('weekly', 'monthly', 'occasional', 'interventions', 'public works')
  
        list(means,ses)
      }

      plot.margins <- function(Y,X){
      listit.model  <- listit(Y=Y, X=X)
      
      num <- length(listit.model$b.star)
      sensitive <- listit.model$b.star
      control1 <- listit.model$b.control[1:num]
      control2 <- listit.model$b.control[(1+num):(2*num)]
      control3 <- listit.model$b.control[(1+2*num):(3*num)]
      vcov <- listit.model$vcov
      vcov1 <- listit.model$vcov.control[1:num,1:num]
      vcov2 <- listit.model$vcov.control[(1+num):(2*num),(1+num):(2*num)]
      vcov3 <- listit.model$vcov.control[(1+2*num):(3*num),(1+2*num):(3*num)]
      
      weekly <<- cbind(1,as.matrix(subset(X, X[,1] == 1)))
      weekly.not <<- cbind(1,as.matrix(subset(X, X[,1] == 0)))
      monthly <<- cbind(1,as.matrix(subset(X, X[,2] == 1)))
      monthly.not <<- cbind(1,as.matrix(subset(X, X[,2] == 0)))
      occasional <<- cbind(1,as.matrix(subset(X, X[,3] == 1)))
      occasional.not <<- cbind(1,as.matrix(subset(X, X[,3] == 0)))
      interventions <<- cbind(1,as.matrix(subset(X, X[,4] == 1)))
      interventions.not <<- cbind(1,as.matrix(subset(X, X[,4] == 0)))
      public.works <<- cbind(1,as.matrix(subset(X, X[,5] == 1)))
      public.works.not <<- cbind(1,as.matrix(subset(X, X[,5] == 0)))
      
      vals = find.coefs(sensitive, control1, control2, control3, vcov, vcov1, vcov2, vcov3)
      
      par(family="Times")
      coefplot(t(vals[[1]])[,1], t(vals[[2]])[,1], CI=2, ylim=c(-.7, .7), xlim = c(1,5.5), vertical=FALSE, pch.pts = 15, 
               var.las = 1, cex.pts = 1.5, cex.var =1, mar = c(1,5.1,5.1,1), ylab="Change in predicted probability", main = " ",
               varnames = c('            weekly\n            patrols', '            monthly\n            patrols',
                            '            occasional\n            patrols', '        interventions',
                            '            public\n            works'))
      coefplot(t(vals[[1]])[,2], t(vals[[2]])[,2], vertical=FALSE, pch.pts = 16, cex.pts = 1.5, cex.var =1, add=TRUE, offset = 0.1)
      coefplot(t(vals[[1]])[,3], t(vals[[2]])[,3], vertical=FALSE, pch.pts = 17, cex.pts = 1.5, cex.var =1, add=TRUE, offset = 0.2)
      coefplot(t(vals[[1]])[,4], t(vals[[2]])[,4], vertical=FALSE, pch.pts = 18, cex.pts = 1.5, cex.var =1, add=TRUE, offset = 0.3)
      legend('top',  legend = c("use trial by ordeal", "call police", "call UNMIL", "do nothing"),  bty = "o", 
             ncol = 4, cex = 0.7, pt.cex = 1, pch = c(15, 16, 17, 18))
      }


      
# Replicate Figure 3

      list = read.dta("survey-r3_leaders.dta", convert.factors=TRUE,convert.underscore=TRUE)
      
      Y.busthap = cbind(list$r3.l.litrbusthap.count, list$r3.l.licobustcophap.num, list$r3.l.licobustunhap.num, list$r3.l.licobustleavehap.num)
      
      R3.C.UNMIL = cbind(list$r3.l.cunmil.week, list$r3.l.cunmil.month, list$r3.l.cunmil.rare,list$r3.l.cunmilpala, list$r3.l.cunmilmaintain)
      R2.C.CONTROLS = cbind(list$r2.l.croad.dist.rainy,   list$r2.l.cfacilities,  list$r2.l.ccapital.offense, 
                            list$r2.l.ccoll.viol, list$r2.c.canycrime, list$r2.l.cfreq.visitpol)

      colnames(R3.C.UNMIL) = c("r3.l.cunmil.week", "r3.l.cunmil.month", "r3.l.cunmil.rare3","r3.l.cunmilpala", "r3.l.cunmilmaintain")
      colnames(R2.C.CONTROLS) = c("r2.l.croad.dist.rainy", "r2.l.cfacilities", "r2.l.ccapital.offense", 
                                  "r2.l.ccoll.viol", "r2.c.canycrime", "r2.l.cfreq.visitpol")

      X.busthap <- cbind(R3.C.UNMIL, R2.C.CONTROLS)
      
      pdf("Figure3.pdf", 10, 7)
      plot.margins(Y.busthap, X.busthap)
      dev.off()
      
      

# Replicate Figure 4
      
      df_persist = fread("r3_res_freqtype_pref_mlogit_margins.txt", skip = 1, header = T, check.names=T)
      
      rows = df_persist[c(3,5,7,9,11)]$V1
      cols = c("State", "Local", "UNMIL")
      
      r3_c_pref_all = df_persist[c(3,5,7,9,11),
                                 c('r3_c_pref_all','r3_c_pref_all.1','r3_c_pref_all.2')]
      r3_c_pref_all = as.matrix(sapply(r3_c_pref_all, as.numeric))
      rownames(r3_c_pref_all) = rows
      colnames(r3_c_pref_all) = cols
      
      r3_c_pref_all.se = df_persist[c(4,6,8,10,12),
                                    c('r3_c_pref_all','r3_c_pref_all.1','r3_c_pref_all.2')]
      r3_c_pref_all.se = as.matrix(sapply(r3_c_pref_all.se, as.numeric))
      rownames(r3_c_pref_all.se) = rows
      colnames(r3_c_pref_all.se) = cols
      
      pdf("Figure4.pdf", 10, 7)
      par(family="Times")
      coefplot(as.numeric(r3_c_pref_all[,2]), as.numeric(r3_c_pref_all.se[,2]), CI=2, 
               ylab="Change in predicted probability", main = "", mar = c(1,5.1,5.1,1), ylim=c(-.3, .3),
               xlim = c(1,5.5),vertical=FALSE, pch.pts = 15, var.las = 1, cex.pts = 1.5, cex.var =1,
               varnames = c('            weekly\n            patrols', '            monthly\n            patrols',
                            '            occasional\n            patrols', '        interventions',
                            '            public\n            works'))
      coefplot(as.numeric(r3_c_pref_all[,1]), as.numeric(r3_c_pref_all.se[,1]),
               vertical=FALSE, pch.pts = 16, cex.pts = 1.5, cex.var =1,
               add=TRUE, offset = 0.1)
      coefplot(as.numeric(r3_c_pref_all[,3]), as.numeric(r3_c_pref_all.se[,3]),
               vertical=FALSE, pch.pts = 17, cex.pts = 1.5, cex.var =1,
               add=TRUE, offset = 0.2)
      legend('bottom', 
             legend = c("rely on informal", "rely on formal", "rely on UNMIL"),
             bty = "o", ncol = 3, inset = 0.05, cex = 0.9, pt.cex = 1, pch = c(15, 16, 17))
      dev.off() 
      
      
      
      